Density functional calculations of nanoscale conductance 
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Density functional calculations for the electronic conductance of single molecules are now common. 
We examine the methodology from a rigorous point of view, discussing where it can be expected 
to work, and where it should fail. When molecules are weakly coupled to leads, local and gradient- 
corrected approximations fail, as the Kohn-Sham levels are misaligned. In the weak bias regime, XC 
corrections to the current are missed by the standard methodology. For finite bias, a new method- 
ology for performing calculations can be rigorously derived using an extension of time-dependent 
current density functional theory from the Schrodinger equation to a Master equation. 
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I. INTRODUCTION AND NOTATION 

Single molecules used as building blocks such 
as diodes, transistors, or switches have attracted 
much interest basis for a future molecu- 

lar electronics^. Many groups world-wide have 
been performing either experiments or calcula- 
tions. There has been tremendous progress, es- 
pecially in the areas of metallic wires B, [j, 0, 0, 0] 



r— 1 peciaiiy m tne areas oi metallic wires [Ay, |4j, |o|, |g 
}=fl and nanotubes @, I, I, M, El, Hill. How 
ever, comparison between theory and experi- 
ment has been much less successful for molecular 
electronics, i.e., organic molecules between two 
electrodes. Experimentally, obtaining consis- 
tently reproducible results from device to device 
has been problematic [bH 15| . Theoretically, the 
challenge is finding a method to quantitatively 
determine device characteristics with neither 
em piri cal input nor over-parameterization (l6l. 
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In recent years, density functional the- 
ory (DFT) calculations of electronic transport 
through single molecules have been published by 
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an ever increasing number of research groups. 
We focus here on purely electronic transport cal- 
culated with DFT-based methods. The most 
prominent method is the Landauer-type scat- 
tering formalism 2l|, [22|, [23j, |24j, formulated in 
terms of Greens functions in combination with 
ground state DFT. It can be derived using the 
Keldysh non-equilibrium Greens-function for- 
malism [25|, 27|]. In the following, we will call 
this method the standard approach. It can also 
be obtained from elementary scattering theory 
(e.g. [17, 26]), o r usi ng the Kubo linear response 



formula 28, 29, 1131 . However, such derivations 



are limited to at best Hartree-interacting elec- 
trons, as we discuss. Apart from calculating 
the current-voltage characteristics of a coher- 
ent molecular junction in the Landauer scat- 

EL 



tering picture 



26|, |27J, 130|, 131|, 132|, 133|, [35J 



36l . |37I . l38l . [39j, |40j, several additional aspects 
such as electro n-phonon coupli ng 4ll. |52|. [5J 
conformation-induced swit ching [54L ~ 55 , 



581 ] . or interaction with light [55] have also been 
addressed within these methods. Questions dis- 
cussed also include the influence of electrode- 
molecule bond geometry 17. 30. J59I . 60 . l6ll . 62 
or 



effects of a gate electrode [l9l. |42|. 58. 63 



Whereas ground state DFT has become quite 
reliable for calculating the electronic structure 
and other properties of molecules and solids [64]]. 
this success has not extended to transport calcu- 
lations through organic molecules [l7l. 19| . While 
these calculations were originally greeted with 
much enthusiasm, researchers in the many-body 



community have always been skeptical [65j, [66 



671 ] . A simple question to ask is, can such cal- 
culations correctly describe the Coulomb block- 
ade regime? That the answer is patently no has 
raised doubts about the validity of this approach 
among that community. These doubts have been 
further compounded by the fact that calculations 
within this scattering framework usually over- 
estimate the experimentally measured conduc- 
tance of organic molecules by about one order of 



magnitude [17j, 1{|. Only for transport through 
atomic metallic wires[2|, 0, 0, [y], do calculations 
yield results in agreement with experiment, but 
this is not a true test of the method, as the same 
result is found in any calculation yielding a unit 



of conductance per channel. 

The performance of this standard approach is 
the main subject of this review. As we describe 
below, neither of the traditional tests of DFT 
calculations, i.e., direct comparison with exper- 
iment or benchmark testing against more ac- 
curate theoretical methods for smaller systems, 
are generally available for this problem. On the 
one hand, the experimental situation is often not 
well-characterized, while on the other, these are 
transport calculations with systems of up to sev- 
eral hundred atoms. Alternative treatments are 
either prohibitively expensive or of such a mod- 
elistic nature as to not allow meaningful tests. 
Thus, at present, there is no simple way to know 
when the standard approach is accurate or reli- 
able. Instead, we examine carefully DFT treat- 
ments, and show that the standard approach is 
an approximation to a more general approach 
using time-dependent DFT, and from this per- 
spective, its limitations can be deduced. 

Ground- state DFT is based on rigorous theo- 
rems and so, if correctly applied to a problem, 
using a sufficiently accurate approximate func- 
tional, will produce an accurate result. The pur- 
pose of the present article is to ask two simple 
questions: (a) is the present standard approach 
formulation derivable from the basic theorems of 
ground-state DFT, and (b) if so, are our present 
approximations sufficiently accurate for conduc- 
tance calculations? The answers show a variety 
of deficiencies (e.g. inadequacy of the ground- 
state approximation, approximations made by 
using local functionals) in the present theory and 
we do not yet know how important these draw- 
backs are. We don't know how frequently situa- 
tions are encountered in which these limitations 
are quantitatively significant. 

We discuss here three major issues that need 
to be resolved to improve on the present state 
of transport calculations, (i) The first involves 
the accuracy of ground-state functionals. The 
functionals presently used in the Landauer and 
Kohn-Sham approximations might not capture 
enough of the physics to be useful and more im- 
portantly they might not give good qualitative 
or quantitatively accurate results. The worst 
defect we have found is due to lack of deriva- 
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tive discontinuity in LDA and GGA functionals 
which leads to an artificial level broadening and 
can greatly overestimate the conductance 19 



due to incorrect positions of the resonances if the 
molecule is weakly coupled to the leads [68|, 801. 
(ii) The second issue is the missing exchange- 
correlation (XC) contribution in the Landauer 
formula. Present calculations entirely miss this 
contribution to the current. Some groups have 
sought to improve on the issues delineated in (i) 
and (ii) by calculating the XC corrections to the 
current via the gradient expansion corrections 
in the Vignale-Kohn approximation j20| and the 
exact exchange Kohn-Sham potential with the 
Optimized Effective Potential (OEP) Q. The 
exact exchange potential can also be estimated 
with self-interaction corrections (SIC) ((381 
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where the self-interaction errors in LDA DFT 
calculations are subtracted out. Furthermore, in 
the weak bias limit, a careful application of the 
Kubo response equation coupled with the DFT 
formulation, can be used to find the exact an- 
swer (jjj. (iii) The third issue we address is 
an exact theory for finite bias, since transport 
experiments are often conducted under a finite 
voltage drop. Thus, an exact formula couched in 
DFT terms must be derived for these conditions. 

In the last few years, the DFT computa- 
tional transport community has become aware 



of these issues [121, 113, l20j, l68(, and a variety of 
approaches to overcome them have been sug- 
gested. Many are looking to alternative for- 
mulations, such as configuration-interaction (CI) 
in quantum chemistry [70|] or GW in many-body 
physics[7l|, 22], to include effects that are missed 
in present (standard approach) DFT treatments. 
Such calculations are sorely needed, to test the 
DFT formulations against and learn their limita- 
tions. Accurate wavefunction treatments are ex- 
tensively used in both ground-state and TDDFT 
as benchmarks and to provide insight into func- 
tional development [23] . But since such calcula- 
tions are typically far more expensive than DFT 
calculations, and given the chemical complexity 
of the devices that can be built, there remains 
an over-riding need to develop a reliable DFT 
approach. 

Thus a variety of new DFT-based formula- 



tions of the problem are being developed. One 
discussed here includes using a Kohn-Sham ef- 
fective single particle version of a Master equa- 
tion formulation of transport 2.4J which will be 
discussed in section IV A 21 Using TDDFT, an- 
other approach obtains the current by calculat- 
ing the time evolution of a system consisting of 
a molecule coupled to two finite metallic con- 
tacts and turning on a potential ste p, resulting 
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in two different chemical potentials [2f 
third is to use large finite leads, and watch a 
capacitance discharge All three methods es- 
sentially begin from a static distribution, apply 
some change, and allow the system to evolve to 
a steady, but non-equilibrium distribution. For 
non- interacting electrons in the weak bias limit, 
all agree, both with each other and the stan- 
dard approach, but likely disagree in the general 
case of interacting electrons in finite bias. Under 
certain limiting conditions, such as adiabatic ap- 
proximations to TDDFT, and local approxima- 
tions to ground-state DFT, they yield the same 
results. 

Because of the breadth of topics we cover 
in this review, we have collected the notation 
used in various formulas and expressions in Ta- 
ble [J for easy reference. We use atomic units 
(e 2 = h = m = 1) throughout, unless oth- 
erwise stated. So all energies are in Hartree 
(1H = 27.2eV = 627kcal/mol) and distances in 
Bohr-radii (0.529A). 



II. REVIEW OF THE STANDARD APPROACH 

In this section, we first review standard cal- 
culations that utilize a combination of ground- 
state DFT and the Landauer scattering formu- 
lation (the standard approach). We then look 
at approximations commonly employed in the 
course of calculating the conductance with this 
method. 

A. Landauer scattering formulation 

The Landauer scattering formulation can be 
easily understood in terms of a simple schematic 
for a molecular tunneling device, as shown in 
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Left lead 



Molecule 



TABLE I: Notation for formulas 



symbol 



n(r) 

Vext(r) 

Vxc(r) 
v H (r) 
Vtot(r) 
v s (r) 

Xs(r,r',w) 
Xpro P (r, r',uj) 

i 
a 

q 

a = L/R 

i,3 
a, 13 



description 



Q«(r,>,<) 
go 

*C 

Ha/3 
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T(e) 
n(e) 
/ 

eHOMo{£LUMC>) 

a(r,r',ui) 

S T 
S 



electron density as a function of position 
ext. potential due to nuclei and ext. fields 
exchange-correlation potential 
Hartree potential 

fext(r) +vn(r), total electrostatic potential 

KS potential = w ex t(r) + v H (r) + u X c(r) 

KS density-density response function 

proper susceptibility 

(density-density response function) 

many-body density-density response function 

occupied level index 

unoccupied level index 

transition index of transition i — » a 

( 3?*(r)$ a (r) = KS ground state orbitals 

e a — U 

quantum numbers for (left/right) lead 
electrons 

indices for KS orbitals of the device region 
cartesian indices 

energy of electron in lead a = L/R 

with momentum k 

coupling between leads and molecule 

creation (destruction) operator for electron 

with momentum k in lead a 

creation (destruction) operator for electron 

with quantum numbers n on the molecule 

density of states in lead a 

Fermi-Dirac distribution functions in leads 

transition rate to left(right) lead 

surface Green's function for the leads 
hopping matrix describing the coupling 
between leads and molecule, 
its elements are given by Vk a ,n 
self-energy matrices = rg r r^ 
full advanced (retarded, greater, lesser) 
Green's function for the extended molecule 
unperturbed KS Green's function for device 
TDDFT KS wavefunction of the entire system 
TDDFT KS wavefunction projected on the 
leads a = L,R 

TDDFT KS wavefunction projected on the 

central (molecule) region 

block of the TDDFT KS Hamiltonian with 

a,/3 = L(left), R(right), C (center/molecule) 

position of level in a resonant tunneling device 

width of resonance in a resonant tunneling 

device/ coupling to the leads 

transmission coefficient as a function of energy 

spectral density of states 

occupation of level 

Fermi energy 

HOMO (LUMO) level of device 
conductivity (current-current response 
function) 

density matrix for total system 
reduced density matrix 



AV 
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Right lead 




FIG. 1: Schematic of the potential of a resonant tunneling de- 
vice with the LUMO level of the device molecule sandwiched 
between the Fermi levels of the leads which are shifted relative 
to each other by an applied bias voltage AV. When AV = 0, 
the chemical potential is [if- The LUMO and HOMO levels 
have a small width 7, indicating weak coupling to the leads. 

Fig. [TJ In this cartoon, the electrons are 
non-interacting and the system is one dimen- 
sional. The leads are featureless boxes, and the 
'molecule' consists of states localized to the bar- 
rier region. In the cartoon, fip is the chemical 
potential of the entire system, i.e., in its ground 
state and in the absence of a bias. The molecule 
has been drawn so as to be weakly coupled to 
the leads, so that the levels on the molecule are 
only slightly broadened into resonances of width 
7- 

In the Landauer picture, the applied bias raises 
the potential on the left lead by an amount AV. 
There is now an imbalance in the system. If one 
waits a long enough time, eventually many elec- 
trons would flow from the left to the right, and 
re-establish equilibrium, with a common chemi- 
cal potential. This is not the regime we are inter- 
ested in. Instead, on an intermediate time-scale, 
one assumes a steady current is established, that 
is sufficiently small as to have no effect on the 
reservoir levels. 

The current can then be calculated from the 
two-terminal Landauer formula 21] which, for 
this case, is simply 



2 f°° 
/=-/ deAfT(e) 



where 



Af = /L(e) - / R (e) 

= f( e -^-AV)-f{e-^). (2) 

and T(e) is the transmission coefficient at energy 
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e. The factor of two is for the two spin channels. 

This can be easily understood as follows. Con- 
sidered function of energy, only those states 
in the window between [ip and fip + AV can 
carry a net current. Those below are occupied 
on both sides, those above are unoccupied on 
both sides. Each state of energy e in the win- 
dow will transmit an electron with probability 
T(e), yielding that contribution to the net cur- 
rent. The schematic shown in Fig. [1] will yield 
a current-bias curve like that shown in Fig. [2] 
The differential conductance, dl/dAV, will be 




FIG. 2: Schematic current-voltage characteristic of the reso- 
nant tunneling device displayed in Fig. [1] The onset of the 
current occurs around elumo — [If- The step is broadened by 
the coupling 7. 

strongly peaked in the position of the LUMO. 
The result is even simpler in the zero bias limit, 
as dV — > 0. Then, df = 5(fj,F)dV, so that the 
conductance becomes 



G 



dl _ 2T(/i F ) 
dV ~ 



7T 



(3) 



Thus for non-interacting electrons in a fixed 
potential, the Landauer formula is easily under- 
stood and justified. 

B. Interacting Electrons 

The Landauer formula for non-interacting (or 
at most Hartree-interacting) electrons was later 
generalized to interacting electrons by Meir and 



Wingreen [25|], who formulated an algorithm 
for calculating the current using the full non- 
equilibrium Green's functions for the system. 
They employ a second quantized Hamiltonian 
description for the electrons in the leads, the 



interacting region (molecule), and the coupling 
between them. Initially uncoupled, the coupling 
between the leads and the molecule is turned on 
slowly via the Vk a ,n term in Eq. (TjJ: 



H = ^ e kac{ a c ka + H int ({d j n };{d n }) 

k,aeL,R 

+ Yl (VkaAadn + H.C.). (4) 
k,aeL,R 

Here, k refers to the momentum of an electron 
with energy efc a in the left or right lead, labeled 
by a. The creation and annihilation operators 
are denoted by c^c) and d^(d), referring to the 
leads and the molecule, respectively. Then, us- 
ing the continuity equation for the current, the 
Keldysh formalism for the Green's functions and 
allowing the electrons in the device region to in- 
teract while keeping the electrons in the leads 
noninteracting, they find an expression for the 
current when the leads are at different chemical 
potentials: 

I = lj de(tr{[f L (e)Tfj - f R (e)YfJ\{G r i ^ - G^)} 
+ HK-TZ)G&), (5) 

and rjf = 2nJ2 aeL/R _p a (e)V a>i (e)V: 4 (e) where 
i,j indexes the states in the interacting region 
and a indexes the states in the leads. G r , G a , 
and G < refer respectively to the retarded, ad- 
vanced, and lesser Green's functions. Meir and 
Wingreen [2^] derived a simpler formula for the 
case of proportional couplings (rf = aY\ 
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de[f L (e)-Ue)]Im[tr{TG r }}, (6) 



where T = r L r R /(r L + r R ). This however, is a 
strong limitation due to their reliance on sym- 
metric contacts which is never fulfilled for a re- 
alistic system. Only an atomic point contact 
could satisfy this condition since it requires that 
each orbital on the device couples equally to the 
left and right contact. This restriction can be 



removed 1 171. [75| , resulting in a general formula 
for the current in terms of the non-equilibrium, 
self-consistent Green's function. 
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FIG. 3: Landauer approach: schematic representation of a 
benzene- 1,4-di-thiol molecule between two gold contacts. The 
molecule plus gold pyramids (55 atoms each) constitute the 
extended molecule as used in the DFT calculations for the 
standard approach. 



C. The standard approach 

Because of the difficulties involved in solv- 
ing the full many-body problem for the non- 
equilibrium Green's functions exactly for an in- 
teracting system of many electrons, the Green's 
function in Eq. fl§D is usually approximated with 
the ground-state Kohn-Sham effective single- 
particle Green's function (DFT-NEGF). This 
complicates the simple picture of Fig. [T] some- 
what, as the KS potential changes with the ap- 
plied bias. A self-consistent KS potential must 
be found, which will continuously change from 
being raised by AV on the left, to being at its 
equilibrium level on the right. These changes 
will not be confined solely to the molecule, but 
should die off within one or two Fermi wave- 
lengths into the leads. Thus one must define an 
extended molecule as in Fig. [3j which includes 
those parts of the leads where the KS potential 
differs from its non-biased value. 

Also note that in most calculations, the 
molecule is chemically bonded to the leads. Thus 
its levels will be much broader than pictured in 
Fig. [TJ overlapping with one another, and delo- 
calizing into the leads. 

In addition, approximations to the self e nerg y 

!)■ 
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matrix S# and Sj, are made (e.g. [2f 
The coupling between the right lead and the 
device is described by the hopping matrix tr 



(whose elements are just the coupling terms Vk a , n 
in Eq. and similarly for the coupling be- 

tween the left lead and the device. There is no 
direct coupling between the leads as this would 
cause electrons from the left lead to run into the 
right lead until a global equilibrium was reached. 
Then, the self-energy that encapsulates the ef- 
fects of coupling the left contact to the device 
can be written as: 



(7) 



where is the surface Green's function for the 
left lead. An equivalent expression can be de- 
rived for the right contact. The full Green's func- 
tion for the device region, G, can then be written 
in terms of the unperturbed KS Green's func- 
tion, gp, for the extended molecule (molecule 
plus small parts of the leads, Fig. [3]) as 



G 



So 



Ft- 



The coupling matrices T R ^ are given in terms 
of the self-energy matrices as 



T L(R) 



-ifS 



L(R) 



4 > 

J L(R)> 



23£ 



L(R)- 



(9) 



This self-energy only describes hopping onto 
and off the device from the leads, but ne- 
glects the other processes that can occur in the 
leads. The leads are thus assumed to be non- 
interacting. In this situation, a Dyson equation 
for the device region leads to the formula for the 
current for the case of non-interacting electrons 
derived from the more general expression given 
in Eq. ©: 



7T 



rfe[/ L (e)-/ R (e)]tr{G a r R GT L }. 



(10) 

Identifying the trace with a transmission coef- 
ficient, this is identical to Eq. ([I]). This ap- 
proach has been implemented successfully by 

27|, 51, E, 



mam 



32 



grou ps by now (e.g . 18 J2R 
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M, 0, H m, mmSh and cal 

culations have been performed for a large va- 
riety of molecules, coupled to electrodes of dif- 
ferent materials like gold, platinum and silicon. 
Molecules studied vary from H2 to alkyl-chains, 
to molecules built from aromatic subunits, or 
metallo-organic complexes. With sucessful we 
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refer to a calculation in which the described stan- 
dard formalism has been implemented and ap- 
plied in a correct way. Especially early DFT 
based calculations of molecular conductance of- 
ten had some shortcomings in the implementa- 
tion or application, e.g. the coupling to the 
macroscopic leads was not accounted for in an 
appropriate way, which lead to strong deviations 
in the calculated current. Sucessful implementa- 
tions, when applied to the same molecular sys- 
tem and contact geometry with the same func- 
tional performed by different groups yield sim- 
ilar results. However, the calculated conduc- 
tances usually overestimate the experimentally 
measured ones by about one order of magnitude 



17, 19 



By far the best studied molecule is 
benzene-dithiol coupled to g old contacts via the 
sulfur atomsjrjjjp, H H 0, H H 0, EE 
47, 48, l49l l50l. 59l| . Fig J9] shows results from a 
calculation [l9|] using the geometry in Fig. [3] and 
comparing DFT with HF. Other calculations use 
slightly different geometries, like e.g. planar gold 
surfaces and a super cell. Using the latter geome- 
try, calculations [5 11] employing the Master Equa- 
tion approach yield a zero bias conductance of 
~ 0.4Go, in very good agreement with Transi- 
esta calculations [34] 0.4G ) and the results 
in Fig. [9] (0.3 G ), which both employed slightly 
different geometries. 



D. Limitations of the standard approach 

Transport is inherently a time-dependent non- 
equilibrium problem with current flow which is 
not at all within the domain of validity of a 
ground-state DFT description. Its natural de- 
scription is found within time- dependent (cur- 
rent) DFT. Therefore, using ground state DFT 
to calculate the device Green's function incor- 
porates several uncontrolled approximations and 
errors which will be investigated in sections II III 
and [TV|) . The time-dependent (TD) XC po- 
tential can be very different from its ground- 
state counterpart. First, the step from time- 
dependent DFT to ground state DFT misses all 
dynamic effects, i.e., the adiabatic approxima- 
tion of section HVB1 Also, for example, the par- 



titioned system is assumed to be in equilibrium 
and disconnected from the leads at some initial 
time. Slowly, a finite voltage is turned on which 
shifts the chemical potentials of the leads. This 
necessarily drives the system out of equilibrium 
and so the electron distribution in the leads do 
not follow the equilibrium Fermi distributions. 
This effect is probably small, but its presence 
points to the many problems with this approach. 

Other problems with this approach include the 
incorrect placement of energy levels due to the 
missing derivative discontinuity when a local ap- 
proximation to the XC-functional is employed. 
This leads to an error in the location of the res- 
onance peaks (section llllj) . The failure of local 
functionals to reproduce the derivative disconti- 
nuity also produces resonance peaks that are too 
wide in energy which results in a general overes- 
timation of the conductance. 



III. INADEQUACY OF GROUND-STATE 
APPROXIMATIONS 



In this section, we take the standard approach 
prescription at face value, and assume it would 
give the correct conductance if implemented 
with the exact ground-state density functional. 
(In the next sections, we will show that this 
is unlikely to be true in general.) But we ask 
the simple question: using present standard den- 
sity functional approximations (LDA, GGA, hy- 
brids), will we get accurate results? To answer 
this question, we must first review some facts 
that are well-known in the DFT community. 



A. Exact ground-state DFT 

In the top panel of Fig. H] we show the exact 
density n(r) of the He atom, calculated by Um- 
rigar et al. [83j, using quantum Monte Carlo to 
minimize the energy of a highly accurate wave- 
function. In the bottom panel of Fig. [4] we plot 
both the nuclear (i.e. external potential) poten- 
tial, which is —Z/r in atomic units, and the 
exact Kohn-Sham potential f s (r) for this sys- 
tem. Two non-interacting electrons, inserted in 
this potential, reproduce exactly the density of 
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the top panel of Fig. 01 By the Hohenberg- 
Kohn theorem [§3], this potential (if it exists) is 
unique. All modern Kohn-Sham density func- 
tional calculations 8^] are calculations of these 
fictitious non-interacting electrons in a KS po- 
tential. The goal of much research in DFT is to 
provide ever more reliable approximations to the 
exchange-correlation (XC) energy, E KC [n}. Its 
functional derivative with respect to density pro- 
vides the Kohn-Sham potential via 



^ext 



r) + v H (r) 



where v ex t(r) is the original external potential, 
t> H (r) is the Hartree (or classical or electrostatic) 
potential, and f xc [ra](r) = 5E^ c [n]/Sn(r) is the 
XC potential. In traditional ground-state DFT, 




FIG. 4: Top panel - exact radial density for the He atom. 
Bottom panel - external and exact KS potentials for the He 
atom (atomic units). 

the eigenvalues of the KS potential, e,, have no 
formal meaning, except that the eigenvalue of 
the highest occupied molecular orbital (HOMO) 
is exactly the negative of the ionization en- 
ergy (Koopman's theorem), as can be shown by 
studying the asymptotic decay of the density (86|]. 
However, the lowest unoccupied molecular or- 
bital (LUMO) eigenvalue does not in general 



TABLE II: Exact Kohn-Sham energies for the He atom, i.e. 
the orbitals for the KS potential of Fig. 3] 



orbital 


energy [H] 


Is 


-0.90372436 


2s 


-0.15773164 


2p 


-0.12656995 


3s 


-0.06454705 



match the negative of the electron affinity, i.e., 

(12) 



In He, e l8 



-0.903 H, but e 2s 



Clumo = —0.158 if, while the electron affinity 
^4 = 0. Thus the fundamental gap, I — A, is 
not equal to the KS gap, i.e, the HOMO-LUMO 
energy difference, see Table HT1 

What happens then, as electrons are added to 
the system? This question was answered more 
than 20 years ago. Consider a Hydrogen atom, 
far (say lOA) from a featureless metal surface 
(e.g., jellium), as shown in Fig. [5J and ask what 



Molecule 



Lead 





FIG. 5: Density and potential of a molecule-lead system where 
resonance occurs when the chemical potential of the lead lines 
up with €lumo of the molecule. 

happens to the KS potential function of the 
global chemical potential of the system. Since we 
chose the H atom far from the surface, its energy 
levels pick up a tiny width 7, as they are broad- 
ened into resonances. If \i matches the LUMO 
energy, as shown in the cartoon, electrons would 
occupy that level. But as soon as there's even 
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an infinitesimal occupation, the level must move, 
in such a way as to keep the exact Koopman's 
theorem satisfied. Since the density changes at 
most infinitesimally, the only allowed change (in 
the region of the atom) is a constant jump in the 
KS potential, by exactly the amount needed to 
restore Koopman's theorem for the new HOMO. 
This is shown in Fig. [61 




FIG. 6: Exact KS potentials for H and H with an extra in- 
finitesimal electron, illustrating the derivative discontinuity. 



B. DFT approximations 

The success of most DFT calculations is based 
on good approximations to the XC energy it- 
self, as this determines so many properties of 
the system. The most popular approximations, 



such as LDA|85j, GGA[87|, and hybrids of ex- 
act exchange and GGA, such as B3LYP 88| and 
PBE0I891, have well-known successes and fail- 
But they have the following failures in 



ures. 



common, because they are density functionals: 

• Self-interaction error: none are exact for a 
one-electron system, in which 



E x = -U, E=0 



;i electron) (13) 



They all have poorly behaving potentials 
far from the nuclei. The true KS potential 
decays as — 1/r far from a neutral atom, and 
this contribution is from the exchange po- 
tential. In Fig. 0, we have also plotted the 
LDA potential for the He atom. It decays 




f 5 

FIG. 7: Exact and LDA KS potentials for the He atom. 

far too rapidly, and so its orbitals are far 
too shallow. The HOMO is at -0.5704 H, 
while the LUMO is not bound at all. 

• None contain the derivative discontinuity, 
so their potentials do not jump when the 
particle number passes through an integer. 

(Actually, hybrids have about 1/4 exact ex- 
change, but this is not enough to cure these ills.) 

On the other hand, orbital-dependent func- 
tionals cure all these ills (at least, approxi- 
mately). The original and simplest method 
is the self-interaction correction of Perdew and 
Zunger[90|. and is often used with LDA for 
strongly correlated systems. The corrected 
exchange-correlation term is then given by: 



E. 



LDA-SIC 

xc 



72 



(14) 

where n«(r) = |0j(r)| 2 . This functional is ex- 
act for one electron, decays correctly at large 
r, and its potential jumps discontinuously at 
integer particle number. The dashed line in 
Fig. [8] shows the huge improvement in the po- 
tential compared to LDA for the He case. More 
satisfactorily, there are now many implementa- 
tions of exact exchange within DFT, in which 
the orbital-dependent exch ang e is treated as an 
implicit density functional [9 11] . Such optimized 
effective potential (OEP) calculations are often 
prohibitively expensive for large molecules, but 
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FIG. 8: Perdew Zunger self-interaction corrected KS potential 
for He (dashed line). Exact and LDA KS potentials for the 
He atom for comparison. 

are exact for one electron, have a potential that 
decays as — 1/r, and v s (r) jumps discontinuously 
when a fraction of an electron is added. 

Without some form of self-interaction correc- 
tion, continuous density functionals allow elec- 
trons to self-repel, yielding orbitals that are too 
diffuse, in potentials that are too shallow. 

C. Effect on transport 

The missing orbital effects in approximate den- 
sity functionals can have drastic consequences 
for calculations of conductances. For example, 
the missing derivative discontinuity in local ap- 
proximations of DFT affects the magnitude of 
the conductance and misplaces the resonance 
peaks. This effect is strongest when the molecule 
is coupled weakly to the leads. For the exact 
Kohn-Sham potential as discussed above, the po- 
tential is discontinuous, suddenly shifting by a 
constant while the energy remains continuous, as 
the next unoccupied level begins to be infinites- 
imally occupied (see Fig. [6]). The origin of this 
discontinuity is due to the fact that em mo, the 
Kohn-Sham LUMO (lowest unoccupied molecu- 
lar orbital) for the iV-electron system is not the 
same as the Kohn-Sham HOMO (highest occu- 
pied molecular orbital) for the N + 1-electron 
system, as seen in sec. IIII A[ In local approxi- 
mations such as LDA and GGA, the fractionally 



occupied level shifts and moves continuously to 
the position of the HOMO of the N+l system as 
the next unoccupied level begins to be fraction- 
ally populated. The effects for transport calcula- 
tion are broad resonance peaks in transmission 
which lead to finite conductance even at ener- 
gies off resonance, yielding an overall increase in 
conductance that is unphysical. 

This was first illustrated in ref . [l7| , where Ev- 
ers et al. performed a test calculation on a real- 
istic system, finding transmission coefficients for 
benzene- 1,4-di-thiol covalently coupled to two 
gold clusters using DFT with a GGA, and com- 
paring to Hartree Fock (HF) methods. The cal- 
culations were performed only at zero bias and 
the results for the transmission (found from the 
corresponding Green's functions) are shown in 
Fig. [9j The transmission coefficient at the Fermi 



/ 






■ 1 




/ 








/ 








/ 




— ■ Hartree-Fock 




/ 




— DFT (55) 




/ 


\ 








\ 








\ 








\ 








\ 


















s 








s 












1 







FIG. 9: Transmission coefficient over energy for benzene-1,4- 
di-thiol using DFT in the standard approach within a semi- 
local approximation (GGA) (solid line) compared with HF 
results (dashed line). Fermi energy is ~ —5.1 eVfl7l]. 

energy is reduced by a factor of 100 in the HF cal- 
culation. This is likely due to the GGA orbitals 
being too diffuse, due to the self-interaction er- 
ror. 

A more serious problem is due to the incorrect 
positions of unoccupied levels. To probe an un- 
occupied resonance at zero bias, we can apply 
a gate voltage V g to a double barrier resonant- 
tunneling device (DBRTS) perpendicular to the 
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leads, shifting the LUMO down to the Fermi 
level ep (— li at T — 0). In Fig. [H this sim- 
ply reduces all molecular levels by V g . As a level 
passes through e F as a function of the gate volt- 
age, there will be a peak in the conductance. 
When the resonance starts to overlap with e F , 
the exact KS ground-state potential in the region 
of the molecule will differ significantly from the 
off-resonant situation, as it depends on the occu- 
pation, i.e. it will jump discontinuously by the 
derivative discontinuity. This thereby greatly 
changes the transmission characteristics. The 
transmission peaks are not at the position of the 
unoccupied resonances of the ungated situation. 

For any sharp resonance, the transmission co- 
efficient is given by 

r« - (e - e ,^T (7/2)2 ("0 

where e rcs and 7 are the position and width of the 
resonance. In any self-consistent KS treatment 
(including using the exact ground-state func- 
tional), €r and 7 depend on the ground-state 
density, and therefore on the partial occupation, 
< / < 1, of the resonant level. 
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FIG. 10: Double barrier resonant tunneling system with a 
gate electrode. Zero-bias transmission over gate voltage: 
dashed line is self-consistent approximate functional result, 
dotted line is approximate result for coupling 7 — > 0, and solid 
line is exact result. Here e wmo {N) = 0, e H oMo(A r + 1) = 1 
and 70 = O.lpjJ. 

We will now use a simple model to show how 
the use of smooth, approximate density func- 
tionals produces completely erroneous transmis- 
sions (and hence conductances) as a function 
of V^[19f. Defining the KS spectral function 
Aj(e) = 3(trG s (e)) we can write expressions for 



the spectral density of states, n(e) = A s /ir, as 
well as for the transmission T s = 7v4 s /2. This 
yields a simple linear relationship between n(e) 
and the transmission of such a level, n(e) = 
2T s (e)/(77r). The self-consistent occupation / 
of the level is found from integrating over n(e) 
as 

^/><^ + M^} 

(16) 

The transmission can be obtained by inverting 
Eq. (US]): 

T-\e F ) = 1 + tan 2 Mf(e F ) - 1/2)} . (17) 

For simplicity, neglect any dependence of 7 on 
occupation /, i.e., j(f) = 70, as the actual de- 
pendence is expected to be weak and to have 
only little effect on the transmission peaks. The 
transmission can alternatively be expressed in 
terms of the gate voltage. Setting ep = for 
V g = and assuming a shift of the energy levels 
by —V g due to the applied gate voltage (gate ef- 
ficiency^), we can replace €f by V g in Eq. ([TBI) 
and Eq. (fI7]) . thus describing the transmission 
in terms of an applied gate voltage instead: 

T-\V g ) = 1 + tan 2 Mf(V g ) - 1/2)} . (18) 

Any calculation that has a derivative disconti- 
nuity would give the solid line in Fig. [TUJ In this 
example Ae = e HO Mo(^ + 1) - e LUMO (A0 is sev- 
eral eV. The narrow resonance (width 7 = 0.1) 
is positioned at the energy of Chomo (N + l). On 
the other hand, the dashed line is the result for 
a smooth, (semi-)local density functional. As 
the N+l level gets fractionally filled, the res- 
onance moves continuously from e LUMO (AO to 
Chomo (N + 1), resulting in a smearing out of the 
resonance width. It can be seen that the posi- 
tion of the resonance is displaced in this case. 
It is now centered inbetween the LUMO of the 
N electron system and the HOMO of the N+l 
electron system, assuming e rcs = e LUMO + /Ae 
(i.e. a linear dependence of the potential on the 
occupation for the smooth functional). The reso- 
nance peak should be located at the true HOMO 
of the N + l electron system (the solid line). 

Even in the extreme limit of no width of the 
level (7 — * 0), the resonance in LDA remains 
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broad. For weakly coupled leads where, at any 
occupation, 7 <C Ae, the Fermi level is pinned 
to the resonance (e rcs (/) — > e F ) for / 7^ or 
1. This yields e^ = e LUMO + /Ae and using Eq. 
( TT7|) we obtain the dotted line in Fig. [TUJ Thus, 
in a standard calculation using (semi)local func- 
tional, Eq. (TT7T) always produces a broad peak 
whose width is comparable to Ae, thereby over- 
estimating the total conductance of the device 
- even when 7^0. For the case of a linear 
relation as discussed here, this artificial width is 
just Ae/2. In addition, the resonance position is 
incorrect, being Ae/2 too low. 

It is possible to recover the derivative discon- 
tinuity and so avoid these artifacts with meth- 
ods briefly described in section IIIIB1 The self- 
interaction is effectively removed either approx- 
imately from an LDA description via a self- 
interaction correction (LDA-SIC) or through the 
rather expensive, but more rigorous, OEP meth- 
ods. In this method, the self interaction can 
be removed explicitly, but it is computationally 
more feasible to parameterize the self-interaction 
in terms of its atomic components. What re- 
sults is an orbital-dependent functional that in- 
corporates non-local effects, yielding the deriva- 
tive discontinuity. 

Comparison of I-V curves between LDA and 
LDA-SIC was performed by Toher et. a/[68|, [69 



using a tight-binding calculation within the stan- 
dard approach. The differences in the I-V char- 
acteristics were much less apparent for the case 
of strong coupling, confirming that the missing 
derivative discontinuity is most problematic in 
the limit of weak coupling. A cartoon of this 
effect can be seen in Fig. [11] where the cur- 
rent versus applied bias voltage of a device ef- 
fective single-particle energy level is plotted for 
the LDA case (solid line), and the case for an 
artificial step-like energy (dotted line) which em- 
ulates the derivative discontinuity for the exact 
Kohn-Sham potential. The plot on the left hand 
side reflects the situation of weak coupling to the 
leads (7 = 0.2eV), whereas the plot on the right 
hand side reflects the situation of strong coupling 
to the leads ( 7 = 1.2eV). 
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FIG. 11: Current of a single energy level coupled to two metal- 
lic leads as a function of bias. Left figure corresponds to the 
case of weak coupling and right figure corresponds to the case 
of strong coupling. Solid lines are results for LDA and the 
dotted lines are results using self-interaction corrected LDA 
(LDA-SIC). From Ref. 



IV. WEAK BIAS 

In this section, we discuss — in the weak bias 
limit — the errors made due to the standard ap- 
proach scheme. In this limit, we can use Kubo 
response theory to deduce the exact answer, and 
compare with the standard approach. We also 
give an estimate of the corrections in terms of the 
Vignale-Kohn current density functional. Lo- 
cal functionals not only miss the derivative dis- 
continuity but also the XC contributions to the 
electric field response, which leads to an addi- 
tional overestimation of the conductance. In the 
limit of low bias it is possible to describe trans- 
port with the Kubo response formulation. See- 
ing how the adiabatic local density approxima- 
tion (ALDA) fails in this formulation provides 
clues as to the source of the problem and how to 
correct it. A more detailed discussion of calcu- 
lations in this limit is given in section (IV A II) . 

A. Time-dependent DFT 

The usefulness of ground-state DFT has been 
augmented by the development and imple- 
mentation of TDDFTj92|. The Runge-Gross 
theorem [93(1 shows that, under appropriate con- 
ditions, the time- dependent potential t> ext (r,t) 
is a functional of the time-dependent density, 
n(r,t). This allows one to construct time- 
dependent Kohn-Sham equations, and use lin- 
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ear response theory to find TDDFT corrections 
to the KS transitions, making them the true ex- 
citations of the system. 

How the linear response theory works can be 
seen in the density change to an applied per- 
turbation varying as exp(iut), which can be ex- 
pressed exactly in different ways: 

8n{ru) = J d 3 r X (r,r',uj)5v cxt (r,uj) (MB) 

= J d 3 rx P rop(r,r',uj)5vtot(r,uj) (EM) 

= J d 3 r Xs (r,r',uj)5v s (r,u) (DFT) 

where t> tot (ru;) is the sum of the external and 
induced (a.k.a. Hartree) potentials, while 



Sv ext (r, u) + 5v H (r, u) + 5v^ 



(19) 

is the Kohn-Sham potential perturbation includ- 
ing the XC contribution. Different susceptibili- 
ties are used in different contexts: \ is the full 
many-body (MB) response function, giving the 
density change in response to the external per- 
turbation; Xprop is the proper or irreducible sus- 
ceptibility, giving the response to the perturbing 
potential of the total electric field, both exter- 
nal and induced (Hartree), used in electromag- 
netism (EM) [94]. Finally, Xs is the Kohn-Sham 
response function, constructed from KS energies 
and orbitals of the ground-state KS potential: 



Xs(rr u) 



2 E 



UJ — UJ q + iO H 



+ c.c.lu 



-UJ) 



(20) 



and 



{q\fW) = Jd 3 r y>r' $,(r)/(r,rO<Mr') 

(21) 

where q is a double index, representing a tran- 
sition from occupied KS orbital % to unoccu- 
pied KS orbital a, u q = e a — e i; and <& 9 (r) = 
0*(r)0 a (r), where 4>i(r) is a KS orbital. Thus 
Xs is completely given by the ground-state KS 
calculation. For example, it is the response of 
the two non-interacting KS electrons sitting in 
the KS potential of Fig. 2, and u q are the 



TABLE III: Singlet transition energies for He, comparison of 
the true transitions with the Kohn-Sham transitions for the 
exact KS potential [9S|. 



transition 


KS value [eV] a 


true transition [eV] 6 


Is -> 2p 


21.146 


21.221 


Is -> 2s 


20.298 


20.613 


Is — * 3s 


22.834 


22.923 



'The differences between the KS eigenvalues obtained using the 
exact potential [98l |. 
'Accurate non-rclativistic calculations from Rcf. 991. 



differences of the orbital energies listed in Ta- 
ble [TT1 By definition, 5v xc (r,u) causes these 
non-interacting electrons to have the same den- 
sity response as the real electrons. Expand- 
ing this around the original ground-state den- 
sity, and requiring the same density response, 
we find a Dyson-like equation relating the true 



and KS susceptibilities[95[]. This is just the RPA 
equation well-known in other areas, but with the 
Hartree interaction modified to include XC ef- 
fects. One can further translate this problem 
into an eigenvalue problem[96], whose approxi- 
mate solution for transition frequencies is [97]: 



10 



0J, 



+ 2<9|- 



+ /*c|?> (22) 



where 



/ xc M(rr',t- If) = 5v^ c {vt)/5n{v't')\ no . (23) 

is called the XC kernel. Thus, the effect of 
TDDFT is to produce corrections to the KS 
transitions to turn them into the true optical 
transitions of the system. If we had the exact 
f xc (r,r',u) for the He atom density, calculated 
perhaps from a traditional wavefunction calcula- 
tion, and inserted it in the full TDDFT response 
equations, we would get exactly the results of 
Table EU 



B. TDDFT approximations 

Most often, ground-state approximate XC- 
functionals are used, even for the TD poten- 
tial, which is called the adiabatic approxima- 
tion. This often produces excellent excited-state 
properties, and transition frequencies typically 
within about 0.2 eV of the true numbers for 
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molecules lOOl ]. A simple example is the 7r — > 7r* 
transition in benzene in which, in a LDA calcu- 
lation, the KS transition is at about 5 eV, but 
the TDDFT (ALDA) correction correctly shifts 
it to about 7 eVjlOll . 

TDDFT has been implemented in many stan- 
dard quantum chemistry codes, and is run 
routinely to extract electronic excitations of 
molecules |l00l ]. However, as the number of cal- 
culations has grown very rapidly, the limitations 
of the scheme with an adiabatic XC approxima- 
tion are being felt. Even early on [1021 ]. it was 
realized that Rydberg excitations would be miss- 
ing within ALDA or AGGA, because of the poor 
quality potentials of the underlying ground-state 
approximation. The LDA potential of Fig. [7] is 
shallow and short-ranged, and so has no Ryd- 
berg series. Exact exchange or SIC function- 
al take care of this j98|; Fig. [8] shows how 
accurate the LDA-SIC potential is in compari- 
son. Also, double-excitations can be shown to be 



missing within any adiabatic approximation [96 



although frequency- dependent kern els can be 
constructed that restore them [l03| . Charge 
transfer-type excitations also fail[104j. 

For solids, development has been slower, as the 
local and semi-local nature of approximate func- 
tionals means that their effect becomes negligible 
in the thermodynamic limit. This can be seen 
from the fact that the XC kernel in ALDA is in- 
dependent of q, the wave vector in the Fourier 
transform of r — r', as q — > 0, but the Hartree 
term grows as l/q 2 . One cure for thi s pro blem is 
to use TD current DFT (TDCDFT) [l05| , whose 
validity is established in the first part of the RG 
theorem [93|. Related to this fact is the notion 
that no local approximation exists in terms of 
the density, but a gradient expansion in the cur- 
rent was constructed by Vignale and Kohn [l06| 
for linear response. More recently, by compar- 
ison with solutions of the Bethe-Salpeter equa- 
tion, accurate many-body approximations to the 
kernel have been constructed, yield ing excellent 
results for excitonic peaks, etc |l07l |. 

The VK approximation is the only current- 
dependent approximation that is well- 
established, and is often applied to problems 
where non-locality and memory (i.e. beyond 



adiabatic) effects are important. Most impor- 
tantly, it often yields finite corrections where 
ALDA gives nothing, such as to the (0,0) 
comp onen t nee ded for the optical response of 



solids jl05l . Il08| . or the over-polar izab ility of 



long-chain conjugated polymers |73l. I109I ]. But, 
given that it is a simple gradient expansion, its 
quantitative acc uracy in any g iven situation is 



open to question [1091 . IllOl . 11 11 



Constitutive relations 



We can relate the conductivity a with the sus- 
ceptibility x described in section (1IV Al) using 
current continuity 

dn „ , s 

^ = -V-j. (24) 

Since Sn(t) = 5ne tu)t , we have: 

Sn=-— V<5j. (25) 
iu 

Then, using the relations between the vector po- 
tential and the field: 



5E = — — = tooSA, 
at 

and the definitions of the current response: 
Sj a (ru) = / dh'J^Xafsirr'^A^r' 



(26) 



P 



(27) 



where Xa/3 is the current-current response func- 
tion, a te nsor, and a a/ g is the conductivity 
tensor 1121 ] . We get an expression relating the 
response functions for conductivity and suscep- 
tibility: 

1 



®aP — —Xa/3 
10J 



{21 



Also, going back to the definition of the density 
response and using the relations given above in 
equations (I26H27I) . we obtain: 



5n(ruj) = — / d i r^2d a a a p(rr'uj)(-d' (5 v(Yv'uj)) 
U J P 

(29) 
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which leads to the relation between the conduc- 
tivity and susceptibility: 



must be taken in the li miting procedu re to ex- 



(30) 



a/3 



As described in section HV Al in TDDFT, there 
are three equivalent, exact formulas - usually 
studied in reference to the polarizabilities of 
atoms in strong fields - applied to different many 
body descriptions and requiring different inputs, 
to describe the density response to an applied 
electric field. 

There is an analogous response function for the 
current response to an external electric field. As 
with the density response formula, the response, 
given by the conductivity a in this case, is differ- 
ent in each of the exact expressions. The many- 
body, non-local conductivity a(rr'uj) describes 
the response to the external electric field <5E oxt . 
The proper conductivity a prop is in response to 
the total field E tot = E ext + E H , and the single 
particle Kohn-Sham conductivity a s yields the 
response to the full expression for the electric 
field, E ext + E H + E xc , including the unphysical 
XC contributions E xc . 

Siiru) = /dra(r,r»5E cxt (r, W )(MB) (31) 

dr cr prop (r, r', u) 5E tot (r, u) (EM) 

= j dra s (r,r',u)5E s (r,uj). (DFT) 

D. DC transport from Kubo response 

In the limiting case of weak bias, the response 
can be expanded to first order in the electric 
field. In t his, w e follow the logic of Kamenev 
and Kohn |ll3| . To derive the DC transport 
response, a frequency-dependent electric field is 
applied, and the limit u — > is taken while al- 
ways ensuring that vp/u « L, where L is the 
circumference of the ring. As shown by Kamenev 
and Kohn, this reproduces the Landauer formula 
for weak bias for Hartree-interacting electrons. 
Our work can be regarded simple extension 
of this analysis to DFT. Note that extreme care 
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tract the relevant results [114j, [115 

It can be shown that, as u — > 0, the conduc- 
tivity can be rewritten as the transmission coef- 
ficient familiar from the Landauer formulation. 
This limit has to be performed carefully to ob- 
tain the correct current. It has to be assured 
that the excursion length of electrons in the de- 
vice, given by Zp = 7j; is smaller than the region 
of the density response, l p , which in turn has to 
be smaller than the device dimensions L of the 
extended molecule in the DFT calculation (see 
Fig. ED- 

Using the (DFT) response equation for the 
Kohn-Sham susceptibility and the full expression 
for the field we can obtain the correct current. 
We first rewrite the Kohn-Sham non-local con- 
ductance as 

<T s (rr'cj) = fno(r) <5 (3) (r - r') 1 + R{rr'u)\ /(iuj) 

(32) 

where 



R(rr'uj) 



1 P g (r)P;(rQ 

q 



(33) 



defining P,(r) = 0*(r)V0 j (r) - J (r)V0*(r). 
Here <f>i(v) and e« are the KS orbitals and eigen- 
values, uj a = €i — €j, and q = This result 
can be written more compactly in terms of the 
retarded KS Green's function Gg(rr'e) and the 
corresponding KS spectral density 



AJrr'e) 



-Q^(rr'e)/7r, 



(34) 



just as the regular Xs can be. Thus we find, 
exactly, 



1 



def(e){G r s (v r f ,e + Lu) (35) 



+ (G r s )*(rr', e - u)}W A s (rr' e) 



For small uj, only the imaginary component of 
the KS Green's function contributes to R. Ex- 
pansion in powers of uj yields a term linear in lu, 
and an integration by parts yields the DC con- 
ductance entirely in term s of the spectral density 
at the Fermi energy |ll3j ]: 

a s (uj -> 0) = -7rA s (e F , r, r')V ® VA s (e F , r, r'). 

(36) 
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This result is true for the conductance of non- 
interacting electrons in any single-particle poten- 
tial. Next, we specialize to a Id system, to avoid 
complications. Then, Eq. ( l30l) tells us that, as 
uj — > 0, a s becomes independent of position. 
Thus, <t s from Eq. ( |36|) may be evaluated at any 
choice of z and z'. A n eas y choice is z < and 
z' > 0, and one finds 113l 



0) 



T s (e F ) 



71 



(37) 



Since u s (uj — > 0) is just a constant, it can be 
taken outside the integral of Eq. (|32|) . yielding 



5I(u -> 0) 



5/ = ^(^ tot + 5F> 

7T 

T s (e F 



(31 



7Y 



d s r'(5E ext (u) (39) 



+ 5E K (r'uj)+8E xc (r'uj)) (40) 

where 5V^ t = / dz' 5E cxt (tj) + SE H (z r u) is the 
net drop in total electrostatic potential across 
the device, and 



8V 



xc 



dzSE I)XC (z,u^0) (41) 



is the corresponding drop (if any) in the XC- 
potential. 

Thus, Eq. fT4"2"j) would be exact if the exact V xc 
was properly included. But the implementation 
commonly used in the Landauer formulation of 
molecular electronics corresponds, as we will see 
below in IIVE} to only the Hartree response: 

I rfi+SV 

SI = - deT 3 (e,V)(f L (e) - f n (e)) 

1A<v) ' (LANDAUER) (42) 



■SV t 



71 



tot 



Equations (j38|) and (j42J) are identical, ex- 
cept that the standard approach does not in- 
clude the extra exchange-correlation term, 5V XC . 
This derivation has been recently generalized 
to includ e co rrect averaging over the lateral 
directions |ll7i |. 

E. XC correction to current 

The present implementation of the Landauer 
formulation using local functionals includes the 



Hartree piece of the potential and thus correctly 
includes the charging effects, but it is missing 
the XC piece. To see this, consider the XC 
contribution to the voltage given by Equation 
(j4Tp . Since 5E XC = — V "^xc, this implies that 
SV XC = Sv xc (z — > oo) — Sv xc (z — > — oo). But 
far from the barrier, 5p = 0, and so any local or 
semi-local approximation necessitates that Sv xc 
equals zero far from the barrier. Thus 5V XC = 
when working within these approximations [191 ]. 
Thus ALDA and all other local or semi-local ap- 
proximations miss the non-local interactions of 
the exact XC functional. 

Alternatively, integration of the second equa- 
tion in expressions (|32l) would also give the exact 
result since all three formulas are equivalent and 
exact: 



51 



T 



prop 
7i 



tot 



(43) 



but T prop refers to the full proper transmission 
coefficient which cannot be easily calculated for 
a realistic system, and is not the transmission 
through any single-particle potential. 

TDDFT within a local or semi-local approx- 
imation has been shown to produce erroneous 
results when non-local effects become impor- 
tant, such as in the optical response of solids. 
The Vignale-Kohn functional is a non-local func- 
tional in terms of current density that has been 
successfully applied to situations where non- 
locality cannot be i gnored, s uch as in long conju- 



gated polymers [731. 1 1091 . 11181 ] where local approx- 



imations give overestimates on the static polar- 
izabilities. This non-locality also plays a role in 
the non-equilibrium transport problem as seen in 
sectior JTVDl In the regime of weak bias, Koen- 
topp et al. [3] estimate the size of the XC cor- 
rection to the current in the Vignale-Kohn ap- 
proximation. Since ALDA is a local approxima- 
tion, it misses the non-local interactions of the 
exact XC functional. Inclusion of the viscous 
contribution to 5V XC yields a correction to the 
transmission coefficient that reduces its magni- 
tude: 

5V XC /V w -(1 - T(e F ))T(e F )/407r 1 / 2 4 /2 (44) 

A more explicit expression for the XC correc- 
tion can be calculated. Sai et al. 2(| calculate 
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the dynamical response contribution to current 
flow using the Vignale-Kohn correction in TD- 
CDFT. This dynamical contribution is a viscous 
flow component from the XC field that is missing 
in ground state DFT calculations and gives a fi- 
nite correction to the conductance. A current 
density functional theory is necessary because 
functionals that are dependent on the density 
alone don't contain information about the con- 
stant value of the current. The Vignale-Kohn 
construction has an XC field that has both the 
ALDA XC potential and a term dependent on 
the XC stress tensor which in turn is dependent 
on viscoelastic coefficients and velocity fields. 

A dynamical resistance R dyn arises from the 
DC XC field which increases the total resistance 
of the system, thus acting against the external 
field. 

4 r {d z nf 



3e 2 A r 



V- 



-dz 



ri 1 



(45) 



where a and b are points inside the electrodes, A c 
is the cross sectional area, n is the density, and rj 
is the viscosity. A calculation that includes the 
real part of the stress component of the electric 
field yields a correction of 10% 13 . 



V. FINITE BIAS 

Given the limitations of the standard approach 
already discussed, it has been realized that a 
more fundamental derivation of the conductance 
formula is needed, especially one that lends itself 
to a DFT treatment. For example, in DFT, one 
is not allowed to turn off the coupling between 
the molecule and leads, as the Hohenberg-Kohn 
theorem does not apply in empty space, nor does 
the RG theorem allow for time-dependent inter- 
actions between electrons. 

Several suggestions have been made as to 
how to do this, that might appear quite differ- 
ent. Here we discuss and compare just two of 
these: the Master Equation approach, and the 
TDDFT-NEGF approach. The Master Equa- 
tion requires coupling to a dissipative bath, 
such as the phonons, in order to achieve a 
steady current, while the TDDFT-NEGF ap- 
proach achieves a steady current via dephas- 



ing into the continuum. Moreover, the Mas- 
ter Equation allows for periodic boundary condi- 
tions (PBC's) whereas TDDFT-NEGF uses a lo- 
calized system. Finally, because of this, the Mas- 
ter Equation with periodic boundary conditions 
requires TD current DFT, whereas TDDFT- 
NEGF uses the density as the basic variable. 

A. Master Equation 

1. Periodic Boundary Conditions 

The Landauer formulation, and indeed most 
of the literature on transport, uses the con- 
cept of different chemical potentials on the left- 
and right-electrodes, and assumes some steady 
current-carrying state between them. 

Such a situation is not so easy to realize 
within the basic theorems of density functional 
theory, time-dependent or otherwise. Even 
TDDFT req uires starting from some initial 
wavefunction |ll9d . almost always the ground- 
state wavefunction of some system. But the 
ground state of any system of electrons has at 
most one chemical potential, not two. 

Thus useful DFT descriptions begin with a sys- 
tem in its electronic ground-state, and a single 
chemical potential. So far, only the situation 
in which both electrodes are of the same metal 
have been investigated. Furthermore, to avoid 
the difficulties of having infinite potentials far in- 
side the electrodes, a gauge transformation that 
is standard in solid-state physics is applied and 
a solenoidal magnetic field is imposed on a ring 
of material. A vector potential that is linear in 
time, a = Et then gives rise to a uniform electric 
field on the ring, causing a current to flow. 

The same approach is then applied to finite 
bias, and again avoids the need for two different 
chemical potentials. However, a new complica- 
tion arises, as in the presence of the electric field, 
as u — > 0, if L is kept finite, the electrons will 
accelerate indefinitely around the ring, which is 
not the situation we wish to model. Instead, if 
some coupling to the phonons in the system is in- 
troduced, there will be dissipation, and a steady 
state can develop. It is possible to derive an ex- 
tension of TDCDFT that includes dissipation in 
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FIG. 12: Master Equation schematic - periodic boundary con- 
ditions, magnetic field induces electric field on ring. 



a time- dependent Kohn-Sham Master Equation. 
Note that dissipation is unnecessary in the weak 
bias limit of the previous section, as Joule heat- 
ing is second-order in the perturbation 120l . 



2. Master Equation theory 

A master equation approach has been con- 
structed that introduces dissipation via a quan- 
tum mechanical treatment of the Boltzmann 
equation from statistical mechanics 74| . The 
master equation describes the evolution of the 
density of a system coupled to a heat bath 
and has an analogue with the well-studied prob- 
lems o f opt ical interactions of laser fields with 



matter [12 1| - atomic transitions in the presence 
of electromagnetic fields. Among the advantages 
of this approach is the elimination of the artif- 
ical boundary and contact conditions necessary 
in the Landauer formulation, the inclusion of in- 
elastic processes, and applications beyond the 
steady state situation. 
With this approach, there are no reservoirs 



corresponding to the left and right leads at dif- 
ferent chemical potentials and the voltage drop 
across the barrier is an output of the method 
rather than an input as it is for the Landauer 
approach. Instead of the open boundary con- 
ditions employed by the Landauer approach as 
illustrated in figure (ED, periodic boundary con- 
ditions are imposed [74j such that the system 
is a closed circuit with no exchange of electrons 
with semi- infinite reservoirs (see Fig. [T2l . This 
geometry is a neat way to treat an open system, 
avoiding partitioning an infinite system as in the 
standard approach. 

The approach employs a quantum Liouville 
equation for the total Hamiltonian of the sys- 
tem Ht, which contains the device Hamiltonian 
Ho, the phonon bath R, and the electro n-phonon 
coupling potential V. St is the density matrix 
for the total system. Its time evolution is given 
by 

dS T /dt=-i[H T ,S T \. (46) 

After tracing out the bath degrees of freedom, 
what results is the Liouville equation for the re- 
duced density matrix S along with a term that 
encapsulates the dissipation in the system C[S}: 

dS/dt=-i[H,S] + C[S}. (47) 

The dissipation term describes collisions with 
the heat bath and is given by 



C[S] — ^ ^mn(L nm L mn S + SL 

m,n 

2 L mn SL nm ) , 



n.mJ-'mn 



(48 



where L = c! n c m and T mn are the transition prob- 
abilities obtained through Fermi's golden rule. 
To derive explicit expression, the coupling po- 
tential 



-t. 



+ h.c. 



(49) 



m,n,a 



is treated perturbatively to second order. The 
creation and destruction operators for the elec- 
trons/phonons are c^/c and a) /a respectively. In 
the coupling matrix elements, 7^ n ' the indices 
m and n refer to the electrons and a refers to 
the phonons. T mn is then given by: 



D(e r 



i)\ln 



> e r 
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Dt 



mn\ ' l e 



> e r 



(50) 



The electric field is imposed on the system 
through the addition of a time dependent vec- 
tor potential a(t) in the Hamiltonian. Gauge 
transformations are then performed periodically 
to set the vector potential to zero, otherwise 
the Hamiltonian would grow indefinitely leading 
to numerical instability in the implementation. 
One of the problems with previous attempts to 
use a Master Equation formulation within this 
setup is the apparent current continuity viola- 
tion. But it can be shown that current continuity 
is maintained once the dissipative contribution 
to the current is considered (Z4j]. The equation 
of motion for the time dependent density when 
propagating the system under the Master Equa- 
tion is given by 



d(n(r)) 
dt 



\t=o 



■ V (j(r))+Tr(n(r)C[S]S) (51) 



The last term in this equation is the contribution 
due to the dissipative part of the Master Equa- 
tion and so the total current is then given by the 
standard expression using the current operator 
plus a dissipative part due to the propagation of 
the system under the master equation: 

(Mr)) = (j(r)) + (j D (r)) (52) 

The general many particle formulation must be 
simplified to an effective single particle form to 
be of any use in practice. The many body den- 
sity matrix S that satisfies the Liouville equa- 
tion for a Hamiltonian H, must be rewritten in 
an effective single-particle description such that 
the resulting single particle density matrix S s is 
in terms of the eigenstates of H s . These eigen- 
states are the equilibrium single particle eigen- 
states and are related to S s via the expression 
S s = J2 lm fim\l)( m \- The density functional for- 
mulation which maps a system of interacting 
particles into one of noninteracting particles with 
the same density is a natural direction to pro- 
ceed. In analogy with the Hohenberg-Kohn the- 
orem for ground-state DFT and the Runge-Gross 
theorem for TDDFT, it can be proven that for a 
fixed electron-electron interaction, a given C[S] 
and an initial density matrix So, the potential 



is uniquely determined by its time-dependent 
density n(rt). A single particle form of equa- 
tion (HTj) can be recovered by applying pertur- 
bation theory for a weak interaction between 
the non-interacting electrons and the phonons 
in the bath, tracing out the irrelevant degrees 
of freedom, and adding the Hartree potential to 
the single particle Hamiltonian. What results 
is a single particle form of the Master Equation 
with a Kohn-Sham version of C[S] and a single 
particle form of the density matrix expressed in 
the basis of the equilibrium single-particle eigen- 
states indexed by n, m. If the expansion coeffi- 
cients are related to the many body density ma- 
trix S via f nm = trlSc^Cn], then the single par- 
ticle master equation can be written in terms of 
the single particle eigenstates: 



df, 



—i 



^ - ^ ^\ H n pif)fprn fnpHp m (t)} 

P 

(&nm fnm) ^ ^ (J^np 1" mp) fpp 
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In addition, the parameters in the dissipative 
part of the Kohn-Sham Master Equation (see 
equation (1501) ) can be in principal obtained from 
ground-state DFT linear response calculations, 
thereby eliminating the need for any empirical 
parameters. 

3. Master equation results 

The master equation approach is currently un- 
der development, but some i nitial resu lts for test 
systems have been calculated [IH, [ml . The first 
model tested with the master equation approach 
was a simple 1-D double barrier resonant tunnel- 
ing system (DBRTS). Well known results famil- 
iar to the experimental community were derived 
for a DBRTS and show that the effect of inelastic 
collisions, accounted for in the master equation 
formulation, is important in understanding the 
behavior of these devices. 

Among the approximations made for the 
phonons in this model calculation are, that their 
density of states has a parabolic dependence 
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given by D(u) ~ uj 2 . Furthermore, the cou- 
pling between levels m and n is set to a con- 
stant 7 mn = 7 . In the top two panels of Fig. [131 
results for the potential in the absence of a field 
and in the presence of a field with low dissipative 
coupling 7q are given. Results are similar to the 
results from the Landauer formulation except for 
the following important points. 

The voltage drop across the device is an output 
of the calculation rather than an input and there 
can be seen a small voltage drop across the wire 
as well, but this is soon neutralized by screening 
effects. 
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FIG. 13: Potential for the DBRTS. Top panel - no applied 
external field, middle panel - external field for the case of low 
dissipative coupling, 70. The bottom panel shows the total 
current (solid line), the Hamiltonian current (short dashed 
line) and the dissipative current (dashed line). 

The bottom panel of Fig. [T3] displays the total 
current derived. In keeping with the constraints 
of current continuity, the total current (solid 
line) is constant within numerical error due to 
cancellation between the Hamiltonian current 
(short dashed line) and the dissipative current 
(long dashed line). It should also be noted that 
the dissipative current is larger at the contact 
points, indicating that these are the sites of lo- 
cal Joule heating. 

The Master equation approach also predicts 
the hysteresis effects of a DBRTS as demon- 
strated in Fig. [HI At low dissipative coupling 
70, the hysteresis is much more evident than 



at higher values of 70. For stronger dissipative 
coupling, the resonance peak also shifts towards 
lower voltages. 
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FIG. 14: Hysteresis effects for a DBRTS device as demon- 
strated through I-V plots for two different values of dissipative 
coupling 70. Lower 70 is associated with a more pronounced 
hysteresis. 

In Fig. the electron occupation of the level 
is plotted for different values of the applied bias. 
At higher bias, the distribution deviates from 
equilibrium and exhibits a bump in the tail that 
points to charging of the resonant level. This 
charging is the origin of the bistability and the 
hysteresis observed in DBRTs. The bistability 
arises from the non-linearity associated with the 
Hartree potential. The finite current for voltages 
above the resonance peak stems from the inclu- 
sion of dissipation. In the absence of dissipative 
effects the current would go down to zero the 
moment the resonant level becomes fully occu- 
pied and the broadened level no longer overlaps 
with the Fermi level of the lead. When dissi- 
pation is included, electrons can relax into the 
leads leading to a finite current for bias voltages 
above the resonance. The dissipative coupling 
also controls the size of the hysteresis effect. 

4- Chains of gold atoms 

Further calculations within the Master Equa- 
tion approach have studied electronic transport 
through a 3-atom gold wire sandwiched between 
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conditions the Hamiltonian current dominates. 
It was found that a more physical behavior of 
the I-V characteristics is obtained, when only 
the Hamiltonian current is used to model the 
physical current. Neglect of the dissipative cur- 
rent in this system is further supported by the 
fact that the current continuity violation of the 
Hamiltonian current is weak (see Fig. fET 
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FIG. 15: Electron occupation of the resonant level in a 
DBRTS as given by the diagonal elements of the density ma- 
trix in the steady state. A higher bias leads to nonlinear 
effects as evidenced by stronger deviations from the unper- 
turbed Fermi-Dirac distribution. 
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FIG. 16: Three atom gold chain: supercell geometry show- 
ing the atomic wire connected to two gold electrodes (Aulll 
surfaces). The dark atoms indicate the region where dissi- 
pative coupling is present. Periodic boundary conditions are 
applied in all directions. The lateral interaction between a 
wire and the nearest periodic images has negligible effect on 
the current. 



two Au(lll) surfaces |l23 |. Fig. [16] shows the 
unit cell of the periodic system used in the cal- 
culations. Four layers of gold atoms per side are 
included as the contacts. Dissipation is applied 
in the three outermost atomic layers only. Again, 
the phonon density is assumed to be parabolic 
and the coupling takes the form of 7^ = ■y{i\V\j) 
where 7 sets the strength of the dissipation. 

In this application, where a very small period- 
ically separated cell was used, a very large dissi- 
pative coupling was necessary to force a steady 
state (this was obtained by imposing one quan- 
tum of conductance at once specific value of the 
applied bias) . This leads to an unphysically large 
dissipative current. The strength of the dis- 
sipative coupling can be reduced with increas- 
ing system size, eventually reaching its physical 
value for cells that are large compared to the 
electron-phonon mean free path. Under those 



FIG. 17: Calculated I-V characteristics of a 3-atom gold wire. 
The black dots are the Master Equation results. The dashed 
line indicates the characteristics corresponding to one quan- 
tum of conductance. The error bars reflect the fluctuations of 
the Hamiltonian current in the supercell which measure devi- 
ation from current continuity in the numerical calculation. 

The calculated IV-characteristics, shown in 
Fig. H7] reproduces well the linear behavior mea- 
sured experimentallyjl, H, [1]. In Fig [18] we plot 
the total (external plus induced) potential across 
the device. The potential drops occurs mainly 
across the length of the Au wire. Only a small 
portion of the potential drop occurs inside the 
leads. We would expect no potential drop at all 
inside a perfect metal, the small observed drop 
is due to the dissipation in the three outermost 
atomic layers of the leads. 



B. TDDFT-NEGF 

Another method that avoids the use of two 
external chemical potentials with an artificial 
partitioning and manages to obtain a steady 
current in transport calculations is the exact 
non-equilibrium Green's function approach us- 
ing time-dependent density functional theory 
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be expressed in terms of the lesser Green's func- 
tion: 

n{r,t) = -2iG < (r,t;r,t) (56) 

which in turn can be calculated from equations of 
motion for the non-interacting system. To solve 
the non-interacting problem, the system is parti- 
tioned into the molecule(C), and left/right leads 
(L,R) as in the standard approach (see Sec. Ill Cj) . 
The KS Hamiltonian can then be written as a 
3x3 block matrix, yielding for the time evolution 



FIG. 18: 3-atom gold chain — potential: Total poten- 
tial (including external potential, and induced Hartree and 
exchange-correlation potential) averaged over planes perpen- 
dicular to the wire. The external potential in the supercell is 
given, for illustrative purposes, in the position gauge (which 
does not satisfy the periodic boundary conditions) . The black 
dots indicate the position of the atomic planes in the slab, 
whereas the red dots indicate the atoms of the wire. The 
total potential is essentially flat in the electrodes. The large 
drop across the wire is due to the contact resistance. 



(TDDFT) [75|, |7J, |77j, \±M UM- The sys- 
tem begins in thermodynamic equilibrium in its 
ground-state and the leads and device are cou- 
pled. A time- dependent perturbation is imposed 
deep inside the leads, such that the potential ex- 
hibits a step somewhere inside the molecule. 

The method applies the NEGF Keldysh for- 
mulation to the time-dependent KS equations, 
i.e., to a set of non-interacting particles yielding 
the correct density. The time evolution of the 
system is then described by the KS equations of 
TDDFT 



i%{r,t) = H s (r,t)%(r,t), 



(54) 



where H s (r,t) is the time-dependent KS- 
Hamiltonian, and \I/ S (r,t) is the KS wave- 
function. Using the continuity equation 
d/dtn(r,t) = — Vj s (r, t) for the KS current 
density j s (r, t) this produces the correct time- 
dependent current through the device: 



I(t) 



<i 3 r-f-n(r, t) 
dt 



(55) 



where the integral is over some cross-section 
through the molecule. In NEGF, the density can 
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(57) 



jected onto left/right lead (L,R) and the 
molecule region(C), respectively. 

For the left and right lead (a = L,R), we can, 
using the lead Greens functions g a (see descrip- 
tion in Sec. Ill Cj) , obtain an explicit solution for 



the projected wavef unctions: 

= iga(t, 0)^ Q (0)+ / dt'g a (t, t')H aC ^ c (t') 

Jo 

(58) 

g a is defined via [idjdt — H aa (t))g a (t,t') = 
S(t — f) with the appropriate boundary condi- 
tions gait^t) = —i and g a (t, tf) = 0. 

Using this, we can rewrite the expression for 
the molecule region as 



■9 r , . 



H cc {t)^c{t) + 



+t H Ca ga{t,0)V a {0) 
a=L,R 



dt'E(t,t')^ c {t') 

(59) 



where 



E(M') 



H C L(t)gL(t,t')H LC (t') (60) 
+H CR {t)g R {t,t')H RC {t') 



is the self energy accounting for coupling to the 
leads as described in Sec. Ill CI Thus Eq. (I59"|) 
yields, in principle exactly, the time-evolution of 
the molecular wavefunction in the presence of a 
current through the leads. It is non-Hermitian, 
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as it describes electrons flowing from left to right. 
The solution for the wavefunction of the cen- 
tral region is obtained by propagating an ini- 
tial state, i.e. the ground-state of the extended 
sytem in equilibrium, which is obtained in a fash- 
ion analogous to the standard approach. For the 
actual propagation, transparent boundary con- 
ditions are imposed at the lead interfac es an d a 
generalized Cayley method is used (see [l25 ] for 
details). 

The feasibility of the scheme has been tested 
on simple Id systems. These calculations 
demonstrate the independence of the steady cur- 
rent on the history, and show a variety of fea- 
tures, such as non-monotonic dependence of cur- 
rent on bias, and larger transient currents than 
steady-state currents. First applications of this 
approach include: (i) the stud y of the role of 
bound-states in transport |l26l |. Here oscilla- 
tions in the density and the TDDFT KS po- 
tential have been observed, the system does not 
evolve towards a s teady state, (ii) The cou- 
pling to nucleii [l27| |. However, these are all non- 
interacting problems, and so have not tested the 
procedure when interaction plays a roll. 

An important issue of the formalism is whether 
or not a steady current can arise in the absence 
of any dissipation. For non-interacting electrons 
(e.g. the electrons in the KS system), it has been 
found a steady current develops if 

1. the single-particle Hamiltonian becomes 
time-independent as t — > oo 

2. the electrodes form a continuum of states, 
i.e., are infinite, 

3. local density of states on the molecule is 
smooth. 

These all appear reasonable conditions. Further- 
more, the steady current is independent of the 
history of the turning-on of the potential step, if 
the t — > oo Hamiltonian is. 

The crucial point for the steady current is 
the second one. In the presence of a contin- 
uum of states, even non-interacting electrons de- 
phase and a steady current develops. This is 
the mechanism for achieving a steady current in 
this approach, and makes dissipation to phonons 



or many-body scattering effects unnecessary. A 
continuum of states requires infinite electrodes, 
but these are then implicitly included in source 
and sink terms in the resulting equations. 



C. Master equation versus TDDFT-NEGF 

Both the Master equation approach and the 
TDDFT-NEGF approach go beyond the stan- 
dard approach. Each begins from a situation 
for which we have a basic theorem proving a 
functional- dependence of the potential on the 
density, and from which the Landauer formula 
can be derived, at least in the case of non- 
interacting electrons. They clearly yield differ- 
ent results in cases where their conditions differ, 
such as when dissipation is strong in the Mas- 
ter equation, but it is as yet unknown if they 
differ when applied to the same situation (and 
if they do, which one is 'correct'). In this sec- 
tion, we compare and contrast the two different 
approaches. 

Basic variable: In the Master equation 
approach, because the system experiences a 
solenoidal magnetic field and periodic boundary 
conditions are used, the basic variable is the cur- 
rent density. In contrast, the TDDFT-NEGF 
approach has been developed using the density 
itself as the basic variable. While this is more 
familiar within DFT, the current allows develop- 
ment of simple approximations such as Vignale- 
Kohn, yielding simple corrections to the conduc- 
tance. 

Boundary conditions: Almost all present cal- 
culations are performed on localized systems em- 
bedded between two electrodes, using localized 
basis sets, and this is also true for TDDFT- 
NEGF. The Master equation approach both re- 
quires and allows use of plane-wave codes, and so 
can make it much easier to adapt present solid- 
state codes for use in transport calculations. 

Fundamental theorems: TDDFT-NEGF is 
based on the Runge-Gross theorem, which ap- 
plies only to finite systems. Yet the leads must 
be infinite to produce the required dephasing. 
This is an inconsistency in the approach whose 
implication is debatable. The Master equation, 
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on the other hand, required proving the Runge- 
Gross theorem for the Master equation instead 
of the time-dependent Schrodinger equation, and 
so requires introducing new functionals. These 
may (or may not) reduce to the standard ones 
of TDCDFT in the limit of weak dissipation. 

Need for dissipation: The TDDFT-NEGF ap- 
proach demonstrates that a steady current can 
arise without explicit dissipation mechanisms, 
once the leads are infinite. The Master equa- 
tion approach may well reduce to the same result 
in the limit in which the ring size is very large 
and the dissipation small, but this has yet to be 
demonstrated. If so, they become, in that limit, 
simply two different procedures for finding the 
same result. If not, one might well be correct for 
molecular conductance, the other not. 

Weak bias: With small dissipation, the Master 
equation produces the same conductance in the 
zero-bias limit as the Kubo response j74[, and 
therefore includes the XC corrections to the Lan- 
dauer formula as seen in sectio nTV Dl A similar 
result c an b e derived from the TDDFT-NEGF 
formula |l24| . although couched in DFT terms. 
Thus the formalisms agree in the limit of weak 
bias and small dissipation, even for interacting 
electrons. For the description of Joule heating ef- 
fects and phonon scattering, which are not easily 
incorporated into the TDDFT-NEGF approach, 
the Master Equation formalism provides a nat- 
ural framework. 

VI. SUMMARY 

In this brief, non-comprehensive review, we 
have critically examined the present state-of-the- 
art of DFT calculations of transport through sin- 
gle molecules. Our findings are: 

• Even the steady state of current flowing 
through a molecule is not included by the 
basic theorems establishing ground-state 
DFT. 

• The commonly used approximation of 
ground-state DFT in the Landauer formula, 
which we dub the standard approach, has 
a variety of limitations, making it inexact, 



even if the exact ground-state functional 
were known and used. 

• Standard density functional approxima- 
tions, such as LDA, GGA, or hybrids, are 
insufficiently accurate to treat molecules 
weakly coupled to leads, and likely produce 
large overestimates of the current. This ef- 
fect might be the origin of the overestimates 
relative to experiment. Orbital-dependent 
functionals, such as exact-exchange or 
self-interaction corrected LDA (LDA-SIC), 
should perform much better. 

• The standard approach is only a Hartree- 
level theory for the conductance, and ne- 
glects non-local XC corrections to the con- 
ductance. This is demonstrable in the case 
of weak bias. Either orbital-dependent or 
current-dependent functionals are needed 
to even estimate these corrections. 

• For finite bias, several approaches have been 
developed that are within time-dependent 
DFT frameworks, thus addressing the prob- 
lem of the inadequate ground-state ap- 
proximation. Two of these have been de- 
scribed in this review, (i) The Master equa- 
tion approach which includes dissipation to 
phonons. (ii) The TDDFT-NEGF formal- 
ism whis does not have the need for dis- 
sipation. Connections are being developed 
between the two, and time will tell which is 
more practical, reliable, or relevant. 

Open questions for both new approaches and 
any other DFT treatments include the following: 

• Do they agree with the Kubo response 
weak-bias limit discussed in Sec. IIVDP 
Both new formalisms do this. 

• At finite bias, which effects are included or 
not in each approach? 

• At finite bias, in what limits do they agree 
or disagree with each other and with the 
standard approach? We already know that 
within ALDA, all give the same answer 
as the standard approach, but expect dif- 
ferences with non-local non-adiabatic func- 
tionals. 
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• Does one always approach a steady state, 
and is it unique? The hysterises seen in 
model calculations with the Master equa- 
tion is an example of more than one steady 
solution. 

• Is there a dependence on how the poten- 
tial is turned on? Recently, bound states 
of the molecule have been shown to lead 
to infinitely oscillating contributions to the 
current in model calculations using the 
TDDFT-NEGF approach, when the turn- 
on is non-adiabatic. Do these survive in an 
interacting system? 

• Are there infinite memory effects in 
the time-dependent Kohn-Sham potential? 
These are logically possible, and might even 
be necessary, to reproduce the physics. 

• Exactly what features of these theories are 
needed to reproduce strongly-correlated ef- 
fects such as Coulomb blockade, and is 



there any chance to model such features? 

On the other hand, the present Landauer-type 
calculations (the standard approach) yield the 
correct steady solution to the more sophisticated 
approaches when the functional is local in time 
and space (see Sec UVDp . and this may be suf- 
ficient for many purposes. Only more demand- 
ing calculations (be they time-dependent DFT, 
non-local static DFT, or CI or GW) and better 
characterized experiments can tell us what is im- 
portant to reliable first-principles predictions of 
the conductance of single molecules. 
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